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ABSTRACT 

We propose a model for Diffusive Shock Acceleration (DSA) in which stochastic magnetic fields in 
the shock precursor are generated through purely fluid mechanisms of a so-called small-scale dynamo. 
This contrasts with previous DSA models that considered magnetic fields amplified through cosmic ray 
streaming instabilities; i.e., either by way of individual particles resonant scattering in the magnetic 
fields, or by macroscopic electric currents associated with large-scale cosmic ray streaming. Instead, 
in our picture, the solenoidal velocity perturbations that are required for the dynamo to work are 
produced through the interactions of the pressure gradient of the cosmic ray precursor and density 
perturbations in the inflowing fluid. Our estimates show that this mechanism provides fast growth 
of magnetic field and is very generic. We argue that for supernovae shocks the mechanism is capable 
of generating upstream magnetic fields that are sufficiently strong for accelerating cosmic rays up to 
around 10 16 eV. No action of any other mechanism is necessary. 

Subject headings: turbulence, MHD, shock waves, cosmic rays, scattering, acceleration of particles 



1. INTRODUCTION 

Cosmic rays (CRs), relativistic charged particles with 
energies 10 8 — 10 22 eV, constitute an essential part of as- 
trophysical systems (see Schlickeiser 2003). In galaxies 
they often provide pressure and energy densities compa- 
rable to those of magnetic fields and thermal gas. In 
very dense regions, such as the cores of molecular clouds 
or accretion disks they are the only source of ionization 
that that must be present to allow interaction of mag- 
netic field and the fluid (see McKee & Ostriker 2007). 
The origin of CRs has been a subject of debate from 
the beginning of research in the field (Ginzburg & Sy- 
rovatsky 1964). By now, it is accepted, however, that 
galactic CRs at least up to the "knee" in the spectrum 
just above 10 15 eV are most likely generated primarily by 
strong supernova shocks. 

The modern study of CR acceleration in shocks involv- 
ing so-called diffusive shock acceleration (DSA) began 
with Krymsky (1977) and Bell (1978) test particle mod- 
els. They noted that when particles whose mean free 
paths exceed the thickness of a viscous shock are scat- 
tered upstream and downstream of the shock on ideal 
scatterers moving with the bulk fluid, they gain energy 
each time they return across the shock front. The re- 
sulting steady-state CR distribution is a power-law in 
momentum, (/ cx p~ q ) with a slope q = 4M 2 /(M 2 — 1), 
where M is the shock flow Mach number. This asymp- 
totes to q = 4 for strong shocks. It soon became ap- 
parent, if there is efficient injection of fresh CRs at the 
shock, that CRs can extract significant energy from the 
shocked flow. Since the CRs diffuse ahead of the shock, 
this naturally leads to a pressure gradient upstream of 
the shock transition that smoothly decelerates and com- 
presses flow into the shock, forming a shock precursor. In 
strongly modified shocks, the total velocity jump through 
the precursor can exceed that for a classical adiabatic 
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fluid shock. Since CRs with higher energies usually have 
longer mean free paths, they travel further into the pre- 
cursor than lower energy CRs. Consequently, they "see" 
a stronger shock transition, and are accelerated more 
efficiently. This leads to an upward-concave spectrum 
reaching a somewhat shallower slope than the test par- 
ticle spectrum at large energies (Malkov & Drury 2001 
and ref. therein). 

On the other hand, the scattering events encountered 
by the CRs are not ideal, but depend on some com- 
plex physics determined by details of the local electro- 
magnetic field, which is, in turn, modified by the CRs. 
In a diffusive propagation approximation one must pay 
close attention to this local physics to determine the spa- 
tial diffusion coefficient, D xx , and the momentum diffu- 
sion coefficient, D pp , that go into the diffusion-convection 
equation for the quasi-isotropic CR distribution function 
/: 
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(e.g., Skirling 1975), with some source term added for 
injection and assuming that f{x,p) depends only on one 
spatial coordinate x and the magnitude of CR momen- 
tum, p. This equation uses the so-called "local" system 
of reference, where the particle momentum is measured 
with respect to the rest frame of the fluid. Similarly, 
magnetic field perturbations that are the driver of par- 
ticle scattering should also be defined locally over the 
scales that are sampled by a particle gyrating in the mag- 
netic field 3 (see, e.g., Yan & Lazarian 2004). 

3 Fortunately, the modern theory of strong MHD turbulence, 
which uses a notion of "critical balance" (Goldreich & Sridhar 
1995), is formulated in the terms of local, rather than global field 
(see Cho & Vishniac 2000, Maron & Goldreich 2001, Beresnyak & 
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As the fluid and the high-energy particles couple 
through the electromagnetic field, the problem reduces in 
practice to the study of the generation of magnetic fields, 
their spatial structure and the scattering of particles in 
those fields. This problem, formulated in the most gen- 
eral case (e.g. with the fully three-dimensional Vlasov's 
equation) is very hard to treat. In part this is due to the 
nonlinear nature of the coupling between particles and 
fields. Following the original test particle model a num- 
ber of approaches to include more complete physics have 
been adopted in the literature. One of the most pop- 
ular approaches to capture the feedback of CRs on the 
electromagnetic fields has been to apply the streaming in- 
stability mechanism, where particles, escaping upstream 
along the magnetic field, confine themselves by amplifi- 
cation of resonant waves (Lagage & Cesarsky 1983). This 
approach assumes the existence of the background mag- 
netic field, which is usually taken for simplicity to be 
directed perpendicular to the shock; that is, along the 
shock normal. This field was traditionally assumed to 
be of the order of the background ISM magnetic field, 
Bq ~ 5/iG. 

However, it is easy to show that upstream magnetic 
fields of around bjiG are too weak to provide an efficient 
acceleration of the cosmic rays with energies as high as 
10 15 GeV, as required for supernova shocks to produce 
CRs to the knee. PeV cosmic rays have long mean free 
paths in such a field and have a high probability of es- 
caping the shock, so they are not subject to further ac- 
celeration. This poses a serious problem for the shock 
acceleration of galactic CRs. 

To overcome the problem one can argue that the mag- 
netic field in the preshock region can be much stronger 
than its interstellar value and that the free energy avail- 
able for the shock is sufficient to generate much larger 
fields (Volk, Drury & McKenzie, 1984). The magnetic 
field generation, if pursued through streaming instabil- 
ity, leads naturally to a highly nonlinear stage of the 
streaming instability where SB ^> B . The original clas- 
sical treatment of the instability is not applicable in that 
limit. What happens in the non-linear regime has been a 
subject of much discussion in recent years (e.g., Lucek & 
Bell 2000, Diamond & Malkov 2007, Blasi & Amato 2008, 
Zirakashvili et al 2008, Riquelme & Spitkovsky 2009). 
The current driven instability proposed by Bell (2004) 
has moved recently to the center of this scientific debate. 
The driving electric current of that instability comes from 
drift (streaming) of the escaping CRs. The compensating 
return current of the background plasma leads to a trans- 
verse force on the background plasma that can amplify 
transverse perturbations in the magnetic field. Numeri- 
cal simulations suggest that the initial field strength can 
grow substantially, leading in the nonlinear form to disor- 
dered fields with coherence lengths that depend on field 
strength (Bell 2004, Zirakashvili et al 2008, Riquelme & 
Spitkovsky 2009). 

CR streaming instability in the presence of a relatively 
organized magnetic field is an attractive starting point 
in understanding how CRs may amplify magnetic fields. 
We argue, on the other hand, that the classic treatment 

Lazarian 2009a, b). How the turbulence may be modified in the 
presence of cosmic rays is still of an open question, however (see, 
e.g., Lazarian & Beresnyak 2006). 



of streaming instability can not deal with 5B 3> Bq case, 
i.e., can not generate large fields directly. This is because 
a large perturbed field violates the key assumptions of 
particle dynamics made in establishing the instability. 
Various attempts have been made to "renormalizc" the 
magnetic field so that the perturbation SB becomes a 
new effective -Bo- For instance, Diamond and Malkov 
(2007) considered diffusion of magnetic energy from res- 
onant gyroscales to larger scales due to compressibility 
effects. Among their assumptions were weak coupling 
and advection-diffusion equations for wave number den- 
sities 4 . Like most treatments of this problem, their equa- 
tions were also one-dimensional in space, which required 
relatively simple structure of the cosmic ray density with 
respect to magnetic field direction. The case of a highly 
tangled three-dimensional magnetic field with SB ^> B , 
which is the case we explore here, can not be treated in 
such way. 

Indeed, in following a magnetic field line one will typ- 
ically see regions with alternating direction of the CR 
density gradient. Therefore, streaming CRs will produce 
both forward and backward going waves. This is unlike 
the classic picture where mostly forward-going waves are 
produced. In a turbulent field those regions of the alter- 
nating CR gradient will be non-trivially determined by 
the topology of the field and the details of CR diffusion 
(see Appendix A). 

Is it possible to generalize the streaming instability ap- 
proach to be applicable in a tangled magnetic field by 
arguing that each particle scattering creates a magnetic 
field perturbation on its own gyroscale? The difficulty 
with this is that in the tangled SB ^> Bo field such small 
incoherent perturbations are going to average out, leav- 
ing the basic, residual effect that the CRs apply pressure 
to the fluid. Can substantial magnetic fields be efficiently 
generated in this generic case? This is the subject of the 
present paper. 

In what follows, we argue that there is an alternative 
process that can provide fast magnetic field generation. 
It provides sufficiently strong magnetic field without ap- 
pealing to either classical streaming instabilities or their 
modifications. The CR pressure gradient is the dynami- 
cal agent that forms the shock precursor and also drives 
field amplification 5 . 

We discuss below the generation of a strong mag- 
netic field as the precursor interacts with the interstellar 
medium. We assume that the CR pressure is a smooth 
function applied to the fluid, while the magnetic field 
is generated by purely fluid nonlinear mechanisms. The 
magnetic field, in turn, plays the role of CR scatterer and 
accelerator. The fluid is stirred within the precursor on 
large, precursor-sized, scales by the combination of fluid 

4 The assumption of weak coupling is broken in the case of 
strongly interacting waves, which necessarily appear in MHD tur- 
bulence (Goldrcich & Sridhar 1995), while the Fokker-Planck as- 
sumption is often broken even in the case of weak coupling, e.g., 
when the interaction is mediated by wavevectors of comparable 
magnitude (Galtier et al 2000, Chandran 2005). 

5 This physics clearly depends on prior development of a strong 
CR precursor to the shock. That requires waiting only until the 
shock is able to accelerate CRs to relativistic energies (e.g., Kang 
et al. 2009). Even microGauss level fields, without strong amplifi- 
cation, can do that quickly. For instance standard diffusive shock 
acceleration using Bohm diffusion with a 1 fiG field in a 10,000 
km/s shock produces GeV CRs on a timcscale of about a year. 
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Fig. 1. — PDF of density of the inflow fluid is approximately 
log-normal due to the ambient turbulence in the ISM. Pictured 
is the PDF of density from simulations with different sonic Mach 
numbers (Beresnyak, Lazarian & Cho, 2005). 

density inhomogeneities and CR pressure. Although we 
consider the CR pressure applied to the inflowing fluid 
homogeneous in the direction along the shock, the accel- 
eration of the fluid element is highly inhomogeneous due 
to the density inhomogeneity of the fluid. This drives 
precursor turbulence. The magnetic energy is generated 
on intermediate scales through a small-scale dynamo. 
While the cosmic ray pressure gradient is smooth, the 
generated magnetic field is tangled and effectively sup- 
press the streaming instabilities, at least in their classical 
formulation. In assuming a smooth CR pressure distri- 
bution, we are effectively averaging the CR action over 
small scales. We leave the problem of the fluctuations of 
CR and magnetic pressure on gyroscale to future study 6 . 

2. PRECURSOR-GENERATED VELOCITY FIELD 

Although in a classic hydrodynamic shock no informa- 
tion can travel upstream and create perturbations, the 
upstream diffusion of CRs leading to precursor formation 
makes this possible, since CRs are much faster than the 
inflow. Because CRs couple only diffusively to the fluid, 
their distribution will generally be smoother and only 
slowly reactive to local changes in the fluid. This creates 
a number of new possibilities, such as the Drury acoustic 
instability (Drury 1984; Dorfi & Drury 1985; Rang et al. 
1992), which is the enhancement of compressible pertur- 
bations by the CR pressure gradient. Such instabilities, 
however, will only operate on perturbations during the 
limited time that a fluid element crosses the precursor. 
This crossing time r c will be determined by the flow pro- 
file as 

f° dx 

Tc= / " TT . (2) 

On the other hand, generically, a fluid element passing 
through the precursor has inhomogeneities of its own. 
Some level of density inhomogeneity is always present in 
astrophysical plasmas, including the ISM (Armstrong et 

6 In strongly CR modified shocks the acceleration of high energy 
particles mostly takes place within the precursor. Furthermore, 
in those shocks the CR pressure is predominantly in high-energy 
particles. As we will see in §3 those particles do not have a gyroscale 
with the conventional meaning. 



al. 1995) and stellar winds, media that typically trans- 
mit SNR shocks. These inhomogeneities, covering an ex- 
tended range of scales, are usually associated with MHD 
turbulence 7 (see, e.g., Elmegreen & Scalo 2004, Lazar- 
ian & Opher 2009). Also there are so-called Small Ion- 
ized and Neutral (SIN) structures, which exhibit fluctu- 
ations at the scale from 100 to 1000 AU that are order of 
magnitude larger in amplitude than those expected from 
the simplistic extension of the turbulent cascade to small 
scales 8 (Heiles, 1997). The stellar winds, acting as media 
through which many young supernova shocks propagate, 
are also inhomogeneous and turbulent. 

Another generic property of the ISM and winds is cool- 
ing (radiative, collisional, etc). Cooling keeps tempera- 
tures fairly low, which makes ISM fluids pliable to com- 
pression. The typical Mach numbers in the warm ISM 
are observed to be between 1 and 10. Such strongly com- 
pressive flows are characterized by significant nonlinear 
perturbations of density 9 . The hot ISM is usually sub- 
sonic with respect to ambient turbulence, which suggests 
smaller magnitude of the perturbations of density. 

As inhomogeneous fluid flows into the precursor with 
speed Mo it is gradually decelerated by the CR pressure 
gradient until it reaches the speed u\ at the dissipative 
shock front (Fig. [5]) . This deceleration could create per- 
turbations of velocity of the order uq — u\, a difference 
between ballistic velocity of the high-density region and 
full deceleration of the low density regions. The resulting 
velocity field could not be purely divergent, but should 
be partially solenoidal due to the coupling between com- 
pressible and solenoidal motions on the outer scale of su- 
personic turbulence 10 . Simulations of strongly compress- 
ible turbulence normally show a 2:1 ratio of solenoidal 
to compressive motions, corresponding to 2:1 degrees of 
freedom in velocity. However, this ratio was observed 
in simulations of statistically isotropic stationary turbu- 
lence, while in our problem we have a preferred direction, 
perpendicular to the shock front and a limited time for 
the turbulence to develop. Taking into account our lack 
of knowledge of the precursor highly compressible turbu- 
lence we parameterize the amount of solenoidal motions 
with a parameter A s such as the RMS of the solenoidal 
part of the velocity u s be defines as 

u s = A s (u - ui), (3) 

where A s could be of the order of unity or smaller. We 
designate the characteristic scale of solenoidal perturba- 
tions as L, Normally, the supersonic ISM exhibits large 
density perturbations on a wide range of scales. This is 
due to the fact that those are created mostly by ISM tur- 

7 Turbulence arises due to the fact that microscopic dissipation 
coefficients such as viscosity and magnetic diffusivity are small and 
the Reynolds numbers are huge, which makes laminar flows in space 
a practical impossibility. 

8 The origins of these inhomogeneities are still debated. 

9 A log-normal distribution of density is predicted and observed 
in simulations, see, e.g., Beresnyak, Lazarian & Cho (2005), for the 
supersonic isothermal magnetized turbulence, see Fig. [T] while the 
subsonic case usually has a normal distribution (which is a special 
case of log-normal) 

10 An easy way to visualize the effect of generation of the 
solenoidal motions is to consider the density inhomogeneities as 
random obstacles that interact with the flow. It is natural, that 
the flow get a solenoidal component as a result of such an interac- 
tion. 
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Fig. 2. — Solenoidal motions, excited by CR precursor (the real 
picture is three-dimensional). In the frame of the shock the preex- 
isting perturbations enter the precursor creating both compressive 
and solenoidal velocity perturbations (the last being depicted). 

bulcnt slow sub-shocks. Although each secondary slow 
sub-shock typically has a moderate Mach number (~ 2), 
a random action of many sub-shocks creates regions of 
high over-density or under-density that can in rare in- 
stances exceed factors at least 10 (see, e.g., Beresnyak, 
Lazarian & Cho 2005). More generally, the density per- 
turbation from individual sub-shocks is of order unity. 
The characteristic scale here is the distance between slow 
sub-shocks, which is smaller than the turbulence outer 
scale, and could be of order of parsecs. From a practical 
point of view, we will take L from one of the typical scales 
of the problem, either the above distance between sub- 
shocks or the thickness of the precursor, or dynamical 
length u s t c , whichever is smaller 11 . 

In this paper we only consider turbulence hydrodynam- 
ically amplified within the precursor. We do not con- 
sider post-shock turbulence, which is a fairly well-known 
and better explored phenomenon that exists even in pure 
hydrodynamics without a CR precursor (see, e.g., Gi- 
acalone & Jokipii, 2007). For the purposes of this paper 
a rather simplistic description of the post-shock magnetic 
fields is sufficient. It is well known that strong fast shocks 
amplify magnetic field parallel to the shock front by, ap- 
proximately, the shock compression ratio ui/u 2 (where 
U2 is the post-shock velocity). Using that approximation 
we will conservatively assume that the post-shock region 
has the same magnetic field structure as the precursor, 
but with magnetic field strength larger by a factor of 

v /273u 1 /w 2 . 

3. THE SMALL-SCALE DYNAMO 

11 Here we assumed that CR precursor is already mature, i.e., 
the CR pressure is a considerable fraction of shock pressure and 
the significant part of CR energy is in high-energy particles, which 
makes precursor relatively thick. We leave the bootstrap mecha- 
nism of the precursor to the future study. 
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Fig. 3. — Magnetic field spectrum (dashed), generated by a small- 
scale dynamo induced by solenoidal velocity motions (solid). L* is 
an equipartition scale of magnetic and kinetic motions, it plays a 
central role in particle scattering. Upper panel: magnetic and ve- 
locity fields from simulations (Cho ct al 2009) with different dashed 
lines corresponding to magnetic spectra at different times. 

Three-dimensional solenoidal flows can amplify mag- 
netic fields through the strctch-fold mechanism. This 
refers not only to turbulent flows, which have motions on 
all scales, but even to large scale three-dimensional lam- 
inar flows (e.g., Zeldovich et al, 1984). In other words, 
amplification of magnetic fluctuations on small scales is 
a very generic feature of highly-conducting fluids. We 
will consider a so-called generic small-scale dynamo (tur- 
bulent) in which magnetic fields are amplified by ini- 
tially weakly magnetized hydrodynamic turbulence with 
an energy-containing (outer) scale of L. "Small-scale" 
means that the scales of magnetic fields we are inter- 
ested in are smaller than L. The problem of the so-called 
mean-field dynamo, when large-scale magnetic fields are 
generated by small-scale motions is not considered here, 
because the mean-field dynamo is fairly slow (see Vish- 
niac & Cho 2001), while inside the shock precursor the 
time for the amplification is rather limited, since all per- 
turbations are quickly advected to the dissipative shock. 
The small-scale dynamo has three principal stages - a 
kinematic stage, when magnetic energy grows exponen- 
tially, a linear stage and a saturation stage (see Cho et 
al, 2009). 

The kinematic stage of the dynamo has been studied 
extensively by analytic (Kazantsev 1968, Kulsrud & An- 
derson 1992) and numerical tools. The rising spectrum 
with magnetic energy Es(k) ~ k 3 ' 2 down to magnetic 
dissipation scales was predicted and later observed in nu- 
merical simulations. For astrophysical applications the 
kinematic dynamo is irrelevant, since its characteristic 
saturation timescale is of the order of an eddy turnover 
time of the smallest eddies, which is a tiny number com- 
pared to outer timescales. For our purposes we can al- 
ways assume that the kinematic dynamo is saturated and 
the dynamo is in the linear stage. 

In the linear stage magnetic energy grows linearly with 
time as 



where e is the energy transfer rate of the turbulence, 
which can be estimated as e = pu^/L, and Ad can be 
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called an efficiency of the small-scale dynamo. A typical 
spectrum of the velocity and the magnetic field in the lin- 
ear stage is presented on Fig. [3] At each particular time 
when the dynamo operates, the magnetic field reaches 
equipartition with the turbulent velocity field on some 
scale L* . This scale grows with time. On scales smaller 
than L* magnetic and velocity perturbations form an 
MHD turbulent cascade with a fairly steep spectrum. 
On scales larger than L* , the magnetic field has a fairly 
shallow spectrum and velocity has a Kolmogorov spec- 
trum. 

The law of linear growth can be understood as fol- 
lows. The main cascade of energy is down-scale, but it 
is converted from a purely velocity cascade to an MHD 
cascade at a scale L* . One can imagine that part of this 
energy cascades up (an inverse cascade) in the form of 
magnetic energy. Let us call this fraction Ad- In princi- 
ple, Ad can depend on scale, i.e., Ad(L*). However, by 
an argument similar to Kolmogorov's, if the inverse cas- 
cade mechanism is purely nonlinear, then, in the middle 
of the inertial interval there is no designated scale and, 
therefore, there is no dimensionless combination involv- 
ing L* . Therefore, the function Ad(L*) has to be con- 
stant 12 . This gives a linear growth of energy. The linear 
growth can also be obtained if we assume that it takes 
several turnover times to reach equipartition on each suc- 
cessive step to larger and larger L* 13 . A linear growth 
rate has been measured in Cho et al. (2009) as being 
close to Ad ~ 0.06 which is the quantity we will use in 
this paper. 

The shallow part of the magnetic spectrum between 
L and L* normally has a slope a between and —1, as 
observed in simulations. We will need these constraints 
later, when we describe a model of particle scattering. 
We particularly favor a model with a = —1/2. This 
model assumes that while magnetic fields on a scale L* 
are generated by random eddies at the same scale L* 
and contain most of the magnetic field energy, the larger 
scale fields come from equipartition of magnetic tension 
on scale I > L*. This can be estimated as SB 2 (I) /I, while 
the averaged magnetic tension comes from a number N = 
(l/L*) 3 of independent random eddies on scale L*. This 
will give scalings SB(l) ~ Z -1 / 4 and 



E B {k)k = 5B 2 (k) ~ k 1/2 . 



(5) 



When L* approaches L, the small-scale dynamo enters 
the saturation stage in which the magnetic field grows 
more slowly than in the previous linear stage. The satu- 
ration value of magnetic energy depends slightly on the 
level of the mean magnetic field (Cho et al. 2009). In our 
case the mean magnetic field can be considered negligible, 
as the typical AlfVen velocity of warm ISM (~ 12km/s) is 
much smaller than the shock speed and associated turbu- 
lent speeds (see § 2). For the purpose of this paper, how- 
ever, we won't need a saturation stage, since we have lim- 
ited time available for amplification, r c (see § 2), which 

12 This assumes locality of the small-scale dynamo. It is indi- 
rectly confirmed by the linear growth observed in simulations. 

13 Schekochihin & Cowley 2007 proposed a model assuming that 
equipartition is reached at approximately one turnover time, which 
gives a linear growth of magnetic energy. This would correspond 
to our model with Ad ~ 1. 



is normally not enough to reach the saturation stage 14 . 

From the linear stage growth we derive quantities 
SB* = SB(L*,xi) and L*(xi) that we will need in the 
next section: 



SB 2 (L*, Xl ) = 87rA d er(xi) 



SB* 



L*{ Xl ) 
L 



1/3 



and 



dx 



ci u ( x )' 

L*{x\) = (2A d u s T{x l )fl 2 L- 1 ' 2 . 

PARTICLE SCATTERING AND SECOND-ORDER 
ACCELERATION 



(6) 
(7) 

(8) 
(9) 



In this section we derive D xx and D pp of the fast par- 
ticles in the tangled magnetic fields created by the small- 
scale dynamo. It should be noted that particle dynamics 
is considered in the rest frame of the fluid (see also §1). 

As it turns out, there are three different regimes of 
particle scattering, depending on the particle energy, E, 
(see Fig. 2). The magnetic spectrum, described in Fig.[3l 
corresponds to the characteristic magnetic field on a par- 
ticular scale SB (I) ~ \J E(k)k, which increases 15 with 
decreasing scale as /-a/ 2 - 1 / 2 = l^ 1 / 4 for a = —1/2 from 
L until L* and decreases with scale as I 1 / 3 for I smaller 
than L*. 

If the particle energy is sufficiently low, the particle 
will, to the first approximation, gyrate around a mean 
field. Otherwise, its trajectory will be stochastic. In par- 
ticular, if E « e8B{L*)L* , the particle will be gyrating 
along the mean field of SB* = 5B(L*) with Larmor ra- 
dius of r g 2 = E/eSB* . We will refer to these particles 
as low-energy and designate low energies as region (1) in 
Fig. [H For low energy particles the scattering frequen- 
cies and acceleration will be determined by a turbulence- 
based formulae with a turbulence outer scale of L* (the 
scale with largest magnetic field). We will consider this 
case in detail in §4.2. 

4.1. High energy particle scattering 

If the energy of the particle is higher than eSB(L*)L* , 
there is no gyration and the particle's trajectory is fairly 
stochastic. This is due to the fact that for this particle, 
the vector magnetic field will be partially averaged out 
on larger scales (so that SBi ~ Z -1 / 4 ). Let us assume 
that such a particle experiences a Bohm scattering and 
has a mean free path of 



r gl = (E/eSB*)—(L*yT^ 

= (E/eSB*) 4/3 (L*)- 1/3 > L*. (10) 

Indeed, (a) - on such a scale, a particle will be deflected 
by an angle of the order of unity; (b) - the deflection 

14 Indeed, L/u s < r c and saturation requires many L/u s , since 

Ad«l. 

5 Provided that spectrum is sufficiently shallow, a > — 1 
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Fig. 4. — Scattering coefficients as they depend on energy. The 
scattering has three principal regimes: (1), low-energy scattering, 
which depends on the properties of small-scale MHD turbulence 
(two cases where fast modes are present (dashed) and absent (dot- 
ted) are shown); (2), strong scattering, where particles are scat- 
tered efficiently by the strong magnetic fields generated by small- 
scale dynamo; (3), high-energy scattering, where particles are only 
weakly scattered. 

from the larger scale field will be smaller (if I2 > r g i, 
then the deflection angle eSB(l2)r g \/E < 1); (c) - the 
deflection from some smaller scale l\ is also smaller (it 
is a random walk with r gl /lx steps), and the deflection 

angle is eSBih)! 1 / 2 ^ 2 / E < 1 16 ). Assuming E = pc 
from here on, the D xx and D pp for such particles will be 
determined by velocity perturbations on scale r g \, i.e., 
u a rJ( 3 L-V3; thuS; 

D xx = (E/eSB*) 4/3 (L*)- 1/3 c, (11) 
D pp ^p 2 ulr-^L-^lco,E^ (12) 
Finally, for very high energy particles, such as 

E » eSB*L 1 ^ L {L' , ) 1 ^r = e SB* L 3 / 4 ^*) 1 / 4 , (13) 

the trajectory will be a random walk with small de- 
flections from scale L and an effective mean free path of 
E 2 je 2 8B\h. The D pp will be energy independent, as it 
will be set by the outer scale velocity perturbations of 
u s . Consequently, 

D xx = E 2 c/e 2 SB 2 L, (14) 
D pp = (eSB L ) 2 Lu 2 s /c 3 . (15) 

4.2. Low energy particle scattering 



The scattering of gyrating particles in a strong mean 
field has been studied for a long time (e.g., Jokipii, 1966), 
and an associated so-called quasi-linear theory (QLT) of 
scattering has been formulated. A particular property of 
MHD turbulence — its strong anisotropy — makes parti- 
cle scattering from solenoidal modes (e.g., Alfven) fairly 
inefficient (Chandran 2000, Yan & Lazarian 2002). So, 
the compressive fast mode, which is rather isotropic, has 
been proposed as a more efficient scatterer in such set- 
tings (Yan & Lazarian 2002, 2004). The fast mode, how- 
ever, can be strongly damped in realistic ISM environ- 
ments, which causes scattering through the fast mode to 
depend on the properties of the gas, such as temperature, 
density and ionization fraction. For the purpose of this 
paper we ignore complicated issues of low energy scatter- 
ing and acceleration and provide two simple, contrasting 
cases. 

First, we can assume that the fast mode is fully 
damped and that only solenoidal modes survive on scales 
of L* and smaller. In this case we can neglect QLT con- 
tributions from the Alfvenic and slow modes as they are 
very small. On the other hand, there arc strong pertur- 
bations of the magnetic field on the outer scale of L*. 
In this case particles, regardless of energy (provided that 
E « eSB(L*)L*) are going to be reflected by magnetic 
bottles on scale L* and D xx will be independent of en- 
ergy and equal to L*c, while D pp will be determined by 
the speed of the bottles. Thus, 

D xx = L*c, (16) 
D pp ^p 2 u 2 s (L*)- 1/3 L-^/c. (17) 

Similar expressions can be derived from TTD reso- 
nance in the presence of large magnetic field perturba- 
tions (Yan & Lazarian 2008). Both of these diffusion 
rates smoothly transition to the solutions from previous 
section. 

Alternatively, we can assume that the fast mode is 
not damped and that the scattering and second-order 
acceleration are due to the fast mode. If we assume 
that the amplitude of the fast mode is approximately 
in equipartition with other modes on the outer scale of 
sub- Alfvenic turbulence, L* , wc will obtain expressions 
that also smoothly transition to higher energy expres- 
sions from the previous section; namely, 

D xx = cr 1 g / 2 2 (L*) 1 / 2 , (18) 
D pp =p 2 u 2 s {L*)^L- 2 '\-JI\ (19) 

(see, e.g. Yan & Lazarian 2004). Here we have used 
so-called acoustic turbulence scaling SB ~ I 1 / 4 for the 
isotropic fast mode, as in Cho & Lazarian (2002). 

It is also possible that in the shock acceleration re- 
gions, where the density of CRs is high, the scatter- 
ing is affected by collective effects where compressions 
of magnetic field induce the gyroresonance instability in 
the fluid of compressed CRs as discussed in Lazarian & 
Beresnyak (2006). We do not provide a discussion of this 
more complex case here. 



Provided that the spectrum is falling (a < 0) 



5. HEATING 
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Turbulent heating provided by solenoidal motions will 
be of the order of the turbulent rate, e = pu^/L. How- 
ever, this assumes that the turbulent cascade dissipates 
all of its energy into thermal particles. This might not be 
true for efficient second-order acceleration of low-energy 
particles. In particular, if the spectrum of CRs is suf- 
ficiently soft (steep), the second order acceleration will 
drain energy from turbulence and put it into CRs as par- 
ticles would tend to diffuse to higher energies. On the 
other hand, if the spectrum is sufficiently hard (shallow), 
energy carried by high energy CRs may be able to drive 
turbulence 17 . In principle, one would like to monitor en- 
ergy obtained or lost by CRs due to second-order accel- 
eration/deceleration and adjust turbulent heating rates 
accordingly. 

As a second note we mention that it is reasonable to 
believe that dissipation from secondary shocks (similar 
to ISM turbulent slow shocks) created in the precursor is 
going to be comparable to e. The above considerations 
leave a certain degree of uncertainty in the amount of 
thermal heating in our model. 

6. DISCUSSION 
6.1. Evolution of the ideas on shock acceleration 

The problem of nonlinear diffusive shock acceleration 
(DSA) is a mature area of research with many publica- 
tions since the original linear theory papers of Krymsky 
(1977), Bell (1978), Axford et al (1977) and Blandford 
& Ostriker (1978). Most of this research was motivated 
by the observed CR power-law distribution, and the ap- 
parent robustness of the basic model. The emphasis was 
largely, but not exclusively, on the predicted properties 
of CRs, rather than the physical details of their scatter- 
ing. So as long as the acceleration process was described 
by rather simple and well-grounded convection-diffusion 
equation (1), the details of scattering and fluid dynamics 
have only had to be "reasonable" in order to obtain appli- 
cable results. On the other hand, since many reasonable 
models of the detailed scattering physics were proposed, 
there were, in effect, many models of DSA. One of the 
popular models to account for the scattering is based on 
a linear (SB < B) or "almost linear" (SB ~ B) stream- 
ing instability analysis. The most common application 
of the latter case is so-called Bohm scattering (m.f.p. <~ 
Larmor radius). While these models result in internally 
consistent particle spectra as well as the one-dimensional 
structure of the flow (see, e.g., Malkov 1998; Berezhko & 
Ellison 1999; Blasi 2002; Kang & Jones 2007), the ques- 
tion of physical justification remains. 

Another approach has been to elicit some robust prop- 
erties of the acceleration process without regard to 
underlying particle-fluid interactions (see a review of 
Malkov & Drury 2001 and references therein), although 
these studies emphasized the uncertainties in such im- 
portant quantities as the predicted high energy cutoff of 
accelerated protons. This understanding is critical to a 
resolution of the origins of galactic cosmic rays, and espe- 
cially the so-called "knee" . Furthermore, the higher en- 
ergy frontier is more generally important in astrophysics, 
as it could set limits to physics of what is happening in 

17 This process is qualitatively depicted in § 2. The precur- 
sor can generate velocity fluctuations due to inhomogeneity of the 
precursor pressure by a variety of mechanisms. 



such objects as AGNs or core-collapsed supernovae. In- 
deed, while low-energy CRs can be accelerated in almost 
any source, the highest energy CRs require a combina- 
tion of large magnetic field and large correlation length 
Bil to be contained in the source. Otherwise, they easily 
escape. Present day neutrino experiments, such as ICE- 
CUBE, have set an ambitious goal to peek into the hearts 
of these objects, and obtain high-energy CR properties 
free of uncertainties associated with models of escape and 
propagation. 

Recently, the DSA problem has been reconsidered 
again with regard to the problem of acceleration with 
more realistic scattering (Malkov & Diamond 2006) 
and the fluid dynamics in the presence of both the 
strong streaming instability and strong compressibility 
(Diamond & Malkov 2007). The latter approach uses 
the advection-diffusion equation with one spatial and 
one spectral dimension and assumes weak wave cou- 
pling. This is unlikely to be enough to describe three- 
dimensional compressive MHD turbulence (see the dis- 
cussion in §1). Also, we point out that a classic stream- 
ing instability with a prescribed direction of the CR 
gradient along magnetic field is unlikely to be useful 
in a turbulent, strongly amplified field. Although our 
model is at this point phenomcnological in its treatment 
of strongly compressible turbulence and uses this com- 
pressibility to estimate solenoidal motions that generate 
magnetic fields, we believe that we provided a more real- 
istic description of the magnetic fields generation in the 
preshocked gas. 

Generation of the magnetic fields in the postshock re- 
gion was considered in many publications with both nu- 
merical and phenomcnological means (see, e.g., Cowsik 
& Sarkar 1980, Giacalone & Jokipii 2007, Sironi & 
Goodman 2007, Inoue et al, 2009). Sironi & Goodman 
(2007) estimated a vortical energy that appears after 
GRB afterglow shocks due to preexisting density inho- 
mogeneities. However it assumed that the magnetic field 
reaches equipartition with vortical energy and did not 
discuss the structure of the magnetic field. We argue 
that it is the shock precursor fields which are important. 
Also, the spectrum and the structure of this field are 
paramount for understanding physically motivated scat- 
tering coefficients. Inoue et al (2009) provided a detailed 
two-dimensional numerical study of post-shock turbu- 
lence and dynamo action. Two-dimensional dynamics, 
however is very different from a three-dimensional one 
(see, e.g., Biskamp 2003). 

6.2. Current limitations of the numerical approach 

The full numerical treatment of the highly- 
compressible flows in question is difficult. One- 
dimensional studies of such flows arc meaningless, as 
they are likely to produce a picture which is completely 
different from three-dimensional dynamics. Indeed, 
one-dimensional flows necessarily generate shocks in 
a finite time, and the dynamics are fully dominated 
by those shocks (see, e.g. Suzuki et al 2007). The 
three-dimensional dynamics is more complicated, with 
compressible motions containing only a fraction of 
energy and weak sub-shocks playing some, but not 
necessarily a dominate role. Direct three-dimensional 
simulations of supersonic turbulence are not only 
computationally expensive, but inherently limited in 



a number of ways. The best known limitation is the 
range of scales (typically a useful range of scales in a 
fairly computationally expensive 1024 3 simulation is 
4-200 in grid units). More relevant to the DSA problem, 
however, is the limitation in sonic Mach number, which 
is around 20 for a 1024 3 simulation and is determined 
by having strongly compressible motions (v ~ c s ) on 
the grid cell scale, or having a significant fraction of 
matter accumulated in clumps of the grid cell size. 
The requirement of the DSA problem, however, is not 
M 3 ~ 20 but rather can be M s ~ 1000 or even higher in 
more energetic sources (AGNs, relativistic jets in GRB 
sources). In this situation three-dimensional fluid-PIC 
codes that will aim to describe interacting fluid and 
CRs will suffer from the same limitations as fluid codes, 
but also will be unable to describe a huge spread in 
energies of the CR spectrum, due to limited statistics of 
particles. 

6.3. The source of solenoidal motions 

In §2 we assumed that the solenoidal velocity is a frac- 
tion of the velocity drop along the precursor with the 
main mechanism being the pre-existing density inhomo- 
geneities and the precursor pressure field. We noted that 
we expect this mechanism to be fairly efficient (i.e. A s 
close to unity) as long as preexisting density perturba- 
tions dp/ p are of the order unity. However, it is inter- 
esting to consider the possibility that the initial density 
perturbations are enhanced by the precursor to this level. 
For instance, once the total pressure in the precursor, 
P, is dominated by CRs, the weak coupling between 
fluid fluctuations and CRs leads to regions in the pre- 
cursor where VP • Vp < 0. This is unstable to Rayleigh- 
Taylor-likc instabilities that have been shown to strongly 
enhance turbulence in CR modified shocks (Ryu et al. 
1993). While these effects have been demonstrated, the 
resulting solenoidal turbulence has not yet been quanti- 
tatively evaluated. Therefore, we cannot conclude that 
A s could be considered a small parameter even if the 
inflowing density perturbations are small. 

Aside from inflowing density inhomogeneities, another 
possible mechanism of generating turbulence is related 
to the inhomogeneities in the precursor CR density, or, 
more generally, a dynamic three-dimensional turbulent 
interaction between fluid and CR's, which allows ex- 
change of energy in both directions. This includes in- 
stabilities discussed above. Qualitatively, this effect also 
works to create additional solenoidal motions. However, 
this is a more complicated phenomenon that will be stud- 
ied elsewhere. 

We expect the current instability, which is considered 
in detail in the next subsection, to generate density and 
solenoidal velocity perturbations as well. However, due 
to the limitations of the numerical approach (see the pre- 
vious subsection) the dynamics of density in the nonlin- 
ear stage of current instability is not fully understood. 

6.4. Our approach and current driven instability 

In response to the realization that magnetic fields 
could be amplified significantly compared to the back- 
ground field, a model based on a current-driven insta- 
bility was proposed and tested numerically (Bell 2004, 
Vladimirov et al. 2006; Zirakashvili et al 2008, Riquelme 



& Spitkovsky 2009). In the linear instability stage the 
magnetic field grows fastest on the characteristic scale, 
determined by the initial field Bq, and the current jd, 



I = 1/kc 



cB 



(20) 



where jd is from high-energy CRs that are "rigid" 
enough to have qBg/pc « k c (Bell 2004). The linear 
growth rate depends only on the current according to the 
relation 

id 

1 = 



(21) 



The nonlinear saturation stage is characterized by 
slower growth on larger scales (Bell 2004, Zirakashvili et 
al 2008). This growth, however, advects along with the 
fluid and, as we discussed in § 2, has to be limited by the 
time for flow to cross the precursor. It is also worth not- 
ing that Bell (2004) and Zirakashvili et al (2008) used 
weakly compressible simulations with initial conditions 
that did not have strong perturbations in cither fluid 
density or CR density, which makes it totally different 
from our approach. We believe, there is a good physical 
reason why precursor turbulence is often strongly com- 
pressible and inhomogeneous (see §2). 

One can compare growth rates of magnetic energy, pro- 
vided by the small-scale dynamo (eq. 8) and current- 
driven instability (eq. 21), assuming that the current- 
driven instability is in the stage with 5B ~ B , but the 
linear growth rate still holds. The result is 



dB 2 cur _ 52 j d L 
dB 2 dyn ' cB 



VAO 



(22) 



The second part of the RHS in eq. (|22| can be esti- 
mated from the characteristic va of the ISM (~ lOkm/s) 
and could be as small as 10~ 9 , if the shock speed is large 
(~ lOOOOkm/s). The first part of the RHS of eq. d22J) 
can be interpreted as the ratio of the field created di- 
rectly by the high-energy particle current on scale L to 
the initial field Bq- The robust estimate of the current 
jd, however, seems pretty elusive. Suppose, following 
Riquelme & Spitkovsky (2009), we assume that the cur- 
rent is produced entirely by escaping particles and that 
there is a fixed ratio r] esc w 0.05 between the flux of CR 
energy emitted by the shock and the flux of energy of the 
incoming fluid pu\ h /2 and also assume that a character- 
istic energy of escaping particles is E esc = 10 15 eV. Then 
one can numerically estimate the above ratio. Taking 
L = lpc and assuming u s « 0.5u s h (since we assume that 
the shock is strongly modified, i.e. u s h = uq ~ uq — u± 
and A s is of the order unity), we get 



dBl 



dB 2 dyn 



= 1.6 x 10" 



10 15 eV 

Ever 



( Vesc \ f L 

V0.05/ 



VAO 



\5pG ) \12km/s/ \A s (uq — m) 



0.5u s 



(23) 



This difference in field growths is due to the fact that 
in our model the full pressure (~ energy density) of the 
CRs is behind the force that winds up magnetic fields, 
while the driving force of the current instability model 
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comes only from those CRs that are able to freely stream 
a substantial distance through the flow. The efficiency 
of our model is based on an assumed large value of A s 
parameter (of the order unity) . Also we do not consider 
the feedback of CRs to the full three-dimensional fluid 
dynamics, assuming instead that CR pressure is homo- 
geneous in two directions along the shock. 

6.5. Future work 

The self-consistent treatment of the flow profile and 
acceleration of particles using expressions from § 4 will 
be presented in a future publication. The scattering co- 
efficients assume that efficient scattering and accelera- 
tion will be provided for particles with energies up to 
E 2 - 3 = eSB*L* 1/A L 3 / 4 (see Fig.0}. This energy can be 
estimated taking u s t L, u s sw 10 4 km/s and L s» lpc as 
£■2-3 ~ 3 x 10 17 eV. This energy corresponds to a mean 
free path of the order of L. However, as the acceleration 
efficiency is smaller by a factor of u s /c (Hillas, 1984) the 
maximum acceleration energy will be around 10 16 eV. 
The higher energy particles will be scattered relatively 
less efficiently and likely to form a steeper spectrum. 
This estimate is tentative and needs to be confirmed or 
corrected as the self-consistent treatment of Eq. (1) and 
the evolving structure of the precursor and will be avail- 
able from the future work. 

Finally, we would like to mention that field growth in 
our model could be amended by the consistent descrip- 
tion of the front-running region of the precursor where 
turbulence is not yet developed. At this point, however, 
it is not clear whether this region will be dominated by 
the classic streaming instability, strong compressible ef- 
fects (and, possibly, generation of secondary fast shocks, 
as the magnetic field is not yet amplified to prevent cre- 
ation of those), or the current-driven instability. 

7. SUMMARY 

In order to explain efficient acceleration of high-energy 
CRs in supernova shocks we appealed to magnetic field 
amplification in the shock precursor that is induced by 
the small scale turbulent dynamo. In our picture the ve- 
locity field necessary for such amplification appears hy- 
drodynamically as a result of the strong CR pressure 
gradient acting on an initially inhomogeneous medium in 
the preshock region. We assumed efficient conversion of 
the precursor velocity drop into solenoidal motions due to 
strong density perturbations of the inflowing fluid. We 
also ignored CR inhomogeneities appearing as a back- 
reaction and their nonlinear feedback. We estimate that 
magnetic fields produced by such amplification are able 
to efficiently scatter and accelerate CRs with energies up 
to 10 16 eV. 
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from NASA grants NNG05GF57G and NNX09AH78G, 
as well as the Minnesota Supercomputing Institute. 
AL acknowledges the NSF grants ATM 0648699, AST 
0808118 and the NASA grant NNX09AH78G. 
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Fig. 5. — A 2D slice of the 3D particle tracing experiment, show- 
ing magnetic field projection on the plane (arrows) and particle 
density using grayscale (dark being higher density). 



8. APPENDIX A 

In order to demonstrate that the stochastic field gen- 
erated in the cosmic ray precursor will be hard to treat 
properly in the classic streaming instability approach, 
we ran particle tracing simulations in a stochastic three- 
dimensional field, generated by a small scale dynamo 
(discussed in §3). 

The electromagnetic field was obtained through direct 
three-dimensional numerical simulation of the incom- 
pressible MHD equations with turbulent driving, along 
with a zero mean field and small fluctuations as initial 
conditions. The simulation was similar to MHD2b0h, 
described in Beresnyak & Lazarian (2009b), with the 
exception that it was run for a relatively short time 
while the small-scale dynamo was in its linear stage (see 
§3). The time was chosen so that the equipartition scale 
L* was approximately in the middle of the logarithmic 
range of scales. The electric field was obtained assuming 
v A /c= 10" 5 . 

The particles were injected on one side of the cube and 
their relativistic equation of motion was solved by a hy- 
brid quality-controlling Runge-Kutta ODE solver. The 
particles were injected from the left and escaped to the 
right. Fig. [5] presents a slice of this three-dimensional 
experiment, where magnetic fields are represented by ar- 
rows and the CR particle density by grayscale. We see, 
that aside from the obvious left-to-right global gradient 
we cannot identify any clear structure of particle density 
along the field. This invalidates in this case the assump- 
tion of the classic streaming instability, which requires a 
regular particle density gradient along the field lines to 
produce a regular field-aligned number density current. 
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